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ABSTRACT 


This paper presents and Justifies a new algorithm for solving linear constant- 
coefficient ordinary differential equations. It also discusses the computational 
complexity of the algorithm and describes its implementation in the FORMAC system. 
It concludes with a comparison between the algorithm and some classical algorithms 
for solving differential equations that have been previously implemented. 
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INTRODUCTION 

The class of problems we are solving are those of the type: find a particular 

solution y to the ordinary differential equation L[D]y - f(x), where L(t) is a real 
P 

polynomial a t** + a , t” ^ +. ..-f a, t + a-t denotes r-fold differentiation 

and f(x) is a linear combination of terms of the form x" exp(ttfl3)x, a and 3 
real and n an integer. Note that we may restrict our attention to f(x) which 
are monomials of the above form* by using linear superposition of solutions. Note 
also that we may obtain the solution of L[D]y ■ x*^ exp otx (cox $x or sin Bx) 
by finding the real or imaginary part of the solution of L[D]y ■ x*^ exp (ct+iB)x. 

The general solution of L[D]y • f(x) may be obtained by adding a particular 
solution of this equation to the general solution of L[D]y = 0. The latter is 
determined by the roots of the polynomial L(t), which in general cannot be found 
exactly. We shall not consider this problem further here. 
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ALGORITHM AND JUSTIFICATION 

To present the algorltha for finding a particular solution of 
L[D]y “ f(x) ■ x” exp(o<‘10)x, we consider separately the cases for 6 ■ 0 and 
6 ft 0. 

Case 1 g - 0. 

Step 1 Set M(t) • (t-a)”^^. 

Step 2 Use the extended Euclidean algorithm to determine polynomials A(t) and 
B(t) such that A(t)L(t)/(t-a)’^ + B(t)M(t) - I, where r^O is the 
multiplicity of a as a root of L(t). 

Step 3 Set y_(x) « exp ox d"**^ A[D+a]x”, where D x^ » j! x^'*’’^/(j+r) ! . 

c — p 

Case 2 6 ft 0. 

Step 1 Set M(t) - (t^-2at+a^+e^)"‘*‘^. 

Step 2 Determine A(t), B(t) such that A(t)L(t)/(t^-2at+a^+a^)^ +B(t)M(t) » I, 

where r Is the multiplicity of a±16 as a root of L(t). 

2 2 —r n 

Step 3 Set ■ exp ox{D +6 ] A[D+aJ x exp(iBx), where 

[D^+a^' ^ x^ exp(igx) • exp(i3x) C U . (x), with 

» J 

C . - j!(-l)’'i’^‘*'^/(r-l)!(2B)’''^^ 

^ » j 

U ,(x) - i (r+j-k-l)!(-216)*' x’^'^V(j-k)!(r+k)! . 

k-0 

We now establish the correctness of the algorithm in case 1. 

Theorem The function YpCx) determined by the algorithm in case 1 is a 
particular solution to L[D]y ■ f (x) ■ x** exp ox. 

Proof It is readily verified that M[Djf(x) ■ 0, where M(t) is as 
defined In step 1. Note also that L(t)/(t-a) and M(t) are relatively prime, 
so that there are polynomials A(t) and B (t) as described in step 2 . Mote as 
well that D *^x^ ■ x^ , where D is as described in step 3. 




3 


It then follows that L[D]yp ■ {LlD]/fD-a]^} exp ox D A[I>H*]x” ■ 

{LfD]/[D-a]’^} exp ax A[D+a]x*', using the exponential shift, ■ 

{LfDl/tD-a’’^} exp ox A[D+alx" - {L[D]/[D-a]’^} A[D] exp ox x” » {l - B’[D]M[Dj}f (x) « f (x). 


since M[D]f(x) ■ 0. Q.E.D. 


The following result Is needed to prove the correctness of the algorithm 
in case 2. 

Proposition 1 For all positive Integers r, non-negative integers j, and 
2 2 r i 

real numbers g, [D +g ] exp(igx)*C .U , (x) « x*^ exp(igx), C and 
j (x) being as defined in step 3. 

Proof Set F . (x) ■ exp(iftx) C . U ,(x). The result Is equivalent to 

r.J r»J r»J 

[D^+g^]’^ F (x) ■ x^ exp(igx), for all positive integers r, non-negative 
r » j 

integers j, and real numbers g. The latter result is established by induction 
on r. 

For r = 1, j arbitrary. 


foV] F, ^ (X) ^ ..i 

i.J (2g)J+i 


j+1 


k=0 


(~2ig) X 

(k+1)! 


j! i^'*’^ exp(igx) 

\ 

k k-1 
(-2ig)‘^ x*^ 

(-2iB)*^'*'^ x*^ 

(26)^'*'^ 

L 

k*l 

(k-1)! 

k! 

t 


- (-216) I “ x-' exp(igx) . 

Hence the result is valid for r ■* 1, j arbitrary. 

Suppose now that the result is valid for r £ v, j arbitrary. Then 
^ * (-i/2gv) J and it follows from the Inductive hypothesis that 

[D^+g^]^ j (x) = x'^ exp(lgx). Thus (D+2ig]'^ ^ ^ (x) » x^. In 
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addition, [D^+B^l F , . (x) ■ (i/28v) C , exp .igx) [D+2iB]D U .(x), and 

Vt1»j V,j V+x»j 

it follows from lengthy, but routine, computations that [D-t-2iB]D U .(x) 

v+i.j 

- {(v+j)!x'''Vj!(v-l)!} + 2lBv j(x). Hence 

- (-i/2Bv) j exp(iBx)[IH2i6r‘^^ 

« (-i/26v) exp(i6x)[I>f2iBl'^ {(v+j) !x^“Vj ! (v-1)!} + 2ivU^ ^(x) 

* exp(iBx) 

establishing the result for r»v+l, J arbitrary. Hence the result holds for 
all r and j. Q.E.D. 

With this result, the proof of the correctness of the algorithm in case 2 
is similar to the proof in case 1, hence we will omit it. 

The next result is used to determine the polynomials A(t) and B(t) 
mentioned in step 2 of the algorithm (both cases). It is readily verified by 
induction. For completeness we state it for a Euclidean Domain, which is a 
generalization of the ring of polynomials over a field. 

Proposition 2 Let L and N be relatively prime elements of a Euclidean 
Domain D, and A^, e D be such that A^L + B^N * I. For s a non-negative 

( 2 ®) 

Integer, define ®s+l recursively by « A^(I+B^N' ], 

B » B Then A L + B ^ - I, all s. 

Svl s s s 
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IMPLEMENTATION 

We have implemented the above algorithm as a program in the FORMAC73 system [2]. 
In the program's notation, the input is the polynomial L and the function F. The 
output is a particular solution YP to the equation L{D]Y • F(X), a FORMAL chain 
containing the individual terms of YP, and an Integer giving the length of the 
chain. The functit^i r is required to be of the form exp oX(cos gX or sin gX) 
and X AS used as the independent variable in L, F, and the solution, THe program 
includes a check on the correctness of the particular solution YP obtained, by 
direct substitution Into the original differential equation. 

Operations on coefficients are performed using rational mode arithmetic. 

"sers who find large rational coefficients inconvenient to work with may wish to 
truncate them or convert them to floating point. Programs for doing this are given 
in [6). 

The code for our program follows below. 
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ORtVEP: PROCEDURE OPT 10NS( MAIN ) REORDER; 

/• */ 

/* THIS TEST DRIVER IS GIVEN TO ILLUSTRATE THE METHOD */ 

/* OF CALLING LINODE. FORMAL PARAMETERS ARE PASSED AS «/ 

/♦ PL/I CHARACTER STRINGS OF LENGTH EIGHT. FOR THIS EXAMPLE ♦/ 

/♦ L AND F ARE READ FROM INPUT FILE SYSIN. ♦/ 

/♦ */ 

/* INPUT DATA FOLLOWS //GO.SYSIN DO * CARO. */ 

/* note: imbedded blanks are ignored; input is free form. */ 

/* DATA FOR AS MANY ODES AS DESIRED MAY BE INCLUDED. «/ 

/* ♦/ 

/* example: ♦/ 

/♦ INOUT DATA FOR L<X)s5*X**2 ♦ 3*X -t */ 

/* AND F(X)= X**2 ♦ «EA*(2*X) ♦ SIN(3*X) IS OF THE FORM.* */ 

/* 'SAXT*? ♦ 3*X - 1* «X**2 * fS**t2*X) ♦ SIN(3«X)« */ 

/♦ */ 

/• THE FOLLOWING OCL IS NECESSARY IN ANY DRIVER THAT USES ♦/ 

/* LINODE AS AN EXTERNAL PROCEDURE. */ 

/♦ */ 


OCL LINODE ENTRY(CHAR(8).CHAR(8 ).CHAR(8).CHAR(B).CHAR(an ; 
OCL (LPIN.PIN) CHAR(72); 

ON ENOFILE (SYSIN) GOTO DONE) 

FORMAC_CPTlONS; OPTSET ( EXPNO) ? 

OPTSET (LINFLENGTHS72) ; 

DO WHILE ( • I • B) ; 

GET LIST (LPIN.FIN); 

PRINT_OUT (LX="LPIN"; FXs^FIN**); 

CALL LINODE (*LX*. 'FX*. 'SOLN*. 'LIST*. ‘N'); 

PRINT^OUT (R = S0LN; S=N; TsCHAINfLIST) ); 

FNO; 

done: put SKIPC2) EDIT 1 • END OF 0.0. E. SOLVER* MA); 

END driver; 
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LINOOE: procedure (L*F fVPtTERMSfNTERMS) REORDER; 

/* 

/• this program computes particular solutions of n-th order */ 

/♦ LINEAR ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM LCOlYaF^X) •/ 

/A MHERE L( X)sA(NlAXAAN ♦ A I N- 1)*XA A( N- 1 ) A ••• A A(1)*X A A(0) A/ 

/A where the A(I1 are national constants and F is a (RATIONAL! LlNEARA/ 

/A COMSINATION OF TERMS OF THE FORM XAAN A «EAA(ALFAAX) A G(X) A/ 

/A where G(X)= SIN(BETAAX) or C(X)= COS(BETAAX). a/ 

/• N IS A NON-NEGATIVE INTE * . A/ 

/A ALFA ANO beta ARE RAT10N*.4_ CONSTANTS* A/ 

/A (E*G* F=X. FsCOS(X)# F*X ♦ AEAA ( AA X ) AS I N( A AX ) ) • A/ 

/A ♦/ 

/• FORMAL PARAMETERS ARE PASSED TO LINODE AS PL/I CHARACTER A/ 

/A STRINGS OF length EIGHT. A/ 

/♦ A/ 

/A DATA GOING in; A/ 

/A FIRST argument; polynomial L(XI. a/ 

/* SECOND argument: FORCING FUNCTION F(X>* */ 

/A RESULTS COMING OUT! A/ 

/A THIRD argument: PARTICULAR SOLUTION TO L(OIY=F(X), A/ 

/A FOURTH argument: FORMAC CHAIN CONTAINING THE INDIV- A/ 

/A IDUAL TERMS OF THE PARTICULAR SOLUTION* A/ 

/A FIFTH argument: INTEGER LENGTH OF THE ABOVE CHAIN. A/ 

/A A/ 

/A LINOOE WILL EXPAND L* A/ 

/A A/ 

/A THE FOLLOWING OCL MUST BE IN ANY PROGRAM THAT USES LlNOOiI A/ 

/A AS AN EXTERNAL PROCEDURE: A/ 

/A DCL LINOOE ENTRY(CHARO) •CHAR(8) «CHAR( 8! *CHAR(3) »CMAR( 8) } ; A/ 

/A A/ 

/A LINOOE CONTAINS THE FOLLOWING INTERNAL PROCEDURES: */ 

/A POLYOIV: WHICH DIVIDES A POLYNOMIAL BY ANOTHER POLY- A/ 

/A NOMIAL ANO RETURNS THE QUOTIENT ANO REMAINDER. A/ 

/A re: which COMPUTES THE REAL PART OF A COMPLEX NUMBER. */ 

/A IM: which COMPUTES THE IMAGINARY PART OF A COMPLEX A/ 

/A NUMBER. A/ 

/A POLYOOP: WHICH applies a POLYNOMIAL ARGUMFNT AS A A/ 

/A differential operator TO A SECOND ARGUMENT. A/ 

/A CU: WHICH COMPUTES (COS(BETAAX)AAlASlN(BETAAXnACAU . A/ 

/A note: FORMAC VARIABLES BETA & R USED IN CU ARE A/ 

/A GLOBAL TO LINODE. a/ 

/A DIVMFAC; WHICH DIVIDES A POLYNOMIAL BY A GIVEN FACTOR A/ 

/A AND RETURNS A QUOTIENT ANO MULTIPLICITY. A/ 

/A A/ 

/A ADDITIONAL COMMENTS ON THE ABOVE MENTIONED INTERNAL A/ 

/A PROCFOUReS ARE GIVEN IN THE INDIVIDUAL PROCEDURES. A/ 

/♦ A/ 

/A DEFINITION OF FORMAC VARIABLES USED IN LINOOE*. */ 

/A ALFA: ALPHA. A/ 

/A beta: BETA. A/ 

/A ap: a(I). a/ 

/A bp: b(I>. a/ 

/A B?: TEMPORARY VARIABLE USED IN COMPUTING A<0) 6 B(0). A/ 

/A B3: temporary variable USED IN COMPUTING A(0» & B(0) A/ 

/♦ WHEN BETA IS PRESENT IN F. A/ 

/A CP: used TO CHECK SOLUTION YP BY SUBSTITUTION. A/ 

/A f: a term of FT. A/ 

/A ft: total FORCING FUNCTION AS GiVEN BY THE USER. A/ 

/A Lp; the polynomial L. a/ 

/A l»r: the polynomial lp/npaar . */ 
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/♦ n: n. 

/♦ NP; X-ALPHA OR <X-ALPHA)**2 ♦ 8ETA**2. */ 

/• NTEPMFT: FORMAC equivalent of 00 INDEX FROM 1 TO NARGSFT* */ 

/* N^ERMS: SAME AS NARCSFT* 

/* null: used to INI'^IALIZE beta for THE VACUOUS CASE* */ 

/• Ni: N4-1, 

/♦ Q1 : TEMPORARY VARIABLE USED IN COMPUTING 6(0). ♦/ 

/F 02: TEMPORARY VARIABLE USED IN COMPUTING A(0) & */ 

/♦ B(0) WHEN BETA IS PRESENT. ♦z' 

/♦ r: R. 

/* terms: chain OF SOLUTIONS FOR TERMS OF FT* */ 

/♦ x: INDEPENDENT VARIABLE OF L AND F* X IS USED AS */ 

/♦ THE INDEPENDENT VARIABLE FOR ALL PROCEDURES IN L INODE* ♦/ 

/* YP: solution for one term of FT* *<' 

/F YT: temporary variable USED TO COMPUTE AP(X«-ALFA), ♦/ 

/♦ YTOTAL: sum of particular solutions for all TERMS OF FT* */ 

/♦ 

/♦ FIXED BIN VARIABLES USED AS DO LOOP INDICES IN LINOOE: ♦/ 

/* IS: USED IN COMPUTING AP AND BP* 

/* NARGSFT: number of ARGUMENTS IN FT MHEN FT IS A SUM */ 

/♦ CF TFRMS. 

/• NTERMFT: ranges from 1 TO NUMBER OF TERMS IN FT. ♦/ 

/• 

/* character VARIABLES USED TO PASS FORMAL PARAMETERS TO ♦/ 

/* AND FROM LINOOE: L. F. YP» TERMS. NTERMS. 

FORMAC_OPTIONS; QPTSETC EXPNO) ; 

f_functicn (C'ji re; im; polydop); 

DCL (L.F.YP. terms. NTERMS) CHAR(8>; 

OCL (NARGSFT. NTERMFT, IS) FIXED BIN) 

ATOMIZE (TERMS); 

LET( LP = *'L«; FT = **F”); 

LET( YTOTAL=0); /* WILL CONTAIN SUM OF SOLUTIONS 
IF L0P(FT)=?4 THEN NARG SF T=NARGS( FT ) I ELSE NARGSFT=i; 

/• NARGSFT MAY BE 1 EVLN IF NARGS(FT)>1 */ 

DO NTEPMFT=1 to NARGSFT; /♦ LOOP FOR NR OF TERMS IN FT */ 

LET( NTERMFT = *«NTERMFT”) ; 

IF NARGSFT=l THEN LET( F = FT ) ; ELSE LE T ( F= ARG ( NTERMF T , F T ) ) ; UHlGlNAL PAGE IS 
LFT( N1=HIGHP0W(F,X> ♦ I ); OF POOR QUALITY 

/* NOW EXAMINE F AND COMPUTE ALFA AND BETA ♦/ 

LFT( ALFA=F/REPLACE(F, «E**(SJ*X). I))! 

IF IOENT( AlF a; I ) 

THEN LET ( ALFA=0 ) ; 

ELSE LET( ALF A=REPLACF ( ALFA, »E**(*J*X). »J ) ) ; 

LET{ 8CTA=F/REPL ACE( F, COS(SJ*X), t, SIN($J*X). 1)); 

IF ident(beta; 1 ) 
then do; 

LET( BETA3SNULL); 

LET I NP = X-ALFA); 

/* NOW DIVIDE LP(X) BY POWERS OF NP; LPR IS QUOTIENT 
AND R IS MULTIPLICITY OF NP IN LP */ 

CALL DIVMFAC (‘LP*. • NP* , *LPR*. *R*); 

/« NOW COMPUTE A(0)=AP £■ B(0)sBP USING SINGLE 
STEP EUCLIDEAN ALGORITHM ♦/ 

CALL POLYOIV ( *LPR* . *NP* , *Ql*, *32*); 

LET< AP=1/B2; UP=-QI/B2>; 

end; 

ELSE do; 

LET( BETA=REPLACE ( BETA, COS(SJWX), SJ, £IN($J*X), *J)); 
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LET< NP«X**2 - 2*ALFA*X ♦ ALFA**2 4> BETA**2); 

/* NOW otvios LP bt powers of NP ♦/ 

CALL OlVMFAC ('LP** 'NP** 'LPR* • 

NOW COMPUTE A(0) S B<0} •/ 

/• USING TWO STEP EUCLIDEAN ALGORITHM */ 

CALL POLTOIV I'LPR** 'NP* • *01* t 'Ba*); 

CALL POLVOIV C»NP** 'Ba** •Q2** 'BS*)? 

LET( AP«-02/n3: BP«( l<fQl*02)/B3M 
FNO; 

/* COMPUTE A(S) ♦/ 

00 1S«0 PY t WHILE (2**IS < INTEGERINin; 

LFT( APsAP • (I ♦ BP*NP**C2*A**IS«) I ) ; 

LET( PPs80*9P); 

end: 

/♦ NOW FORM AO+ALFA) (E*P(-ALFA*X)WFIX) ) ♦/ 

LET( YP=«E** (-ALFAWX 1 * F); 

LET( YTsFVALlAP. X« XfALFA) ) I 
LET( YPsPOLYOOPI YT.YP) ) ; 

/♦ NOW APOLY appropriate ANTI-OIFF OPERATOR ♦/ 

IF IOENTI0FTA:NULL) 

then LET! YP=REPLACE{ XAXAYP. XWASJ* 

FAci Aj-a)*x**( sj-a^R)/FACcsj“a»R) ) ) ; 

ELSE IF -^IDENTIRIO) THEN 

LET( YP=REPLACE<X»X*YP. X J*COS ( BET A* X } « RE (CU( S J- 2) ) # 

XP*SJPSINIBETA*X) • IMICUI SJ-2)> ) ) 

LET! YP= WE<»*( ALF A*X> ♦ VPI; 

/* NOW CHECK RESULT */ 

LET( CP=POLVOOP(LP#YP) ) ; 

IF -iI0ENTCCP;F ) then 

oo; 

PUT SKIP! a I EDIT 

C*L(0)YP -a term of F, CP=L(D)YP PRINTED AS DEBUG AIO«)(A) 

print_out( CP) ; return; 

ENO ; 

LET! QUEUE ( TF. RMS ) =YP ) ; 

LET( YTOTAL=YVOTAL*YP) ; 

end; /♦ of do ntermft=i to nargsft; ♦/ 

LET ("YP^sYTOTAL; ••TERMS**sCHA|NITERMS)) ; 

LET ( NTERMS=«NARGSFT«; •*NTERMS“sNTeRMS) ; 

bOLYDIv: PP0C(A,M.Q.REM) reorder; /♦ REM=A (MOO M) •/ 

/* THIS ROUTINE DIVIDES A POLYNOMIAL A BY A POLYNOMIAL M */ 

/• AND RETUPNS THE QUOTIENT Q AND REMAINDER REM */ 

OCL (A.M.O.REM) CHAR(8); 

/♦ THE FOLLOWING ARE TEMPORARY VARIABLES USED IN POLYOIV: */ 

LOCALIZE ( at:ca;cm;na:nm:q;qt:si: 

LET( QT=0; AT = **A**; NA=H 1 GHPOW ( A T • X ) ; NM=Hl GHPOW ( «M«« * X M ; 

IF IOFNT(NM;0) then oo; LET(*«Q**aAT/*«MM;*»REM" = 0) ; RETURN; FNO; 

FUSE LET (CMsCOEFF("M*».X**NM) ) ; 

DO WHILE (INTEGER(NA) >s I NTEGER( NM ) ) ; 

LET( CAsCOEFFIAT. X«WNA)); 

LET( 0=CA/CM * X**(NA-NM}); 

LET( qt=ot+q; s=q*“m"; atsat-s); 

LET( NA=HIGMPOW( AT »X ) ) I 
END ? 

LET( •♦Q« = QT; «REM'*sAT); 

ENO POLVOIV ; 

re; p^_broc <z) reorder; 

/♦ ThTs routine COMOUTES the real PART OF A COMPLEX NUMBER Z */ 
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F_«eTURN (EVALf z. «i «on ; 

P.ENO re; 

ir: f.proc <z) REORoew; 

/♦ THIS ROUTINE COMPOTES THE IMAGINARY PART OF A COMDEX NOMRCR Z */ 
F^RETORN (RlAtEVALlZ.Rl .Ol-Z)): 

f_eno !m; 

POLVOO*»: F_PROC I0,V) REOROERl 

/• THIS ROUTINE RETURNS THE RESULT OF APPLYING THE POLYNOMIAL */ 

/* FIRST ARGUMENT AS A DIFFERENTIAL OPERATOR TO THE SECOND «/ 

/• ARGUMENT. «/ 

F^RETURN (REPLACEIXFXFU* XFFSJ. OER IV( V « X « S J>2 ) ) ) I 
F^ENO POLYDOP; 

cu; F_pRoc (j) local <c;u) reorder; 

/* This routine computes icosibetafxi^rifsinisetaaxdacfu */ 

/* RETA C R ARE GLOBAL VARIABLES . J IS A FORMAC INTEGER. */ 

DCL K FIXED BIN ; 

LET( UsO); 

00 K=0 TO INTEGERIJI; LET( » ; 

LET! U»U ♦ FAC(RFJ-K-t) / FACU-KI / FACIR^K) 

* (-2*Fl*BETA )**K ♦ X**(R^KU; 

END ; 

LET! C=FACIJ) / FAC(R-l) * (-IJFFR 
♦ »!*F(RtJ) / ( 2*BCT A)«F(R4>J ) ) ; 

F_RETURN ( (COS(BETA*X)^*lFSIN(BETAFX) )*C*U) ; 

F_ENO cu; 

OIVMFAC; PRCC<L»»NP,LPRtR) REORDER; 

/* THIS routine divides A POLYNOMIAL LP REPEATEDLY BY A FACTOR ♦/ 

/• NP AND returns THE QUOTIENT POLY AS LPR ANO THE MULTIPLICITY •/ 

/♦ OF NP IN LP AS R */ 

DCL (LP.NP.LPF.R) CHAR(B); 

/• DE'^INE LOCAL TEMPORARY VARIABLES: •/ 

DCL IR FIXED BIN; 

LOCALIZE (orx;llp;rp;qp ); 
let( rp=^o; QP=*«LP"); 

DO IR=-l BY I VHILE < IDENT(RP:0)) ; 

LET! LLP=CR) ; 

CALL PCLYOIV ('LLP'. ‘NP** ‘OP*. ‘RP* ) *. 

end; 

LET! "LPR^sLLP; ORX**»IR»»; •*R"»ORX> ; 

END OIVMFAC ; 


END LI NODE 
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TEST PROBLEMS 


The progreiB was able to solve problems from s variety of sources, 
of typical results follows below. 


A sample 
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5 4 3 2 

-?x ♦ISX -3ftX ♦51X-152 


2 X ;; X 

FX*X #F ♦X C0S(3X)«E tSlNCSXl 


42X 3*X 2 X 

9 * 1/476 X 4E - p/2197 X ffE - 117/13690 X SIN ( 3 X ) 4E 


2 

♦ 1/448 X SIN 1 3 X 

) - 1/13690 

2 X 2 

X COS 1 3 X 1 ♦ 1/312 X COS 

2 

f 3 X ) 64/28461 X 

2 X X 

•E ♦ 4241 1/126E325 X SIN ( 3 X ) «E 

13142/1246X24 X COS ( 3 

X 

X ) *E 4 

2 X 

336/371293 X »E - 21294621/ 


X X 

448540250 SIS ( 3 X ) *F ♦ 9418897/468540250 COS C X X ) *£ 


S s X 


4?X 32X 22X 

T X CM41N ( 1/474 X 4E - 8/2197 X *E 46/28461 X 


♦ X34'371297 X 

2 X 2 X 

*- « - 117/13690 X SIN ( 3 X ) 4F - 1/13690 

2 

X 

COS { X X 

X 

) 

X 

* 42411/1244325 X SIN ( 3 X » 4E - 1 X 1 5 2/ X ’46395 

X 

rns ( ’ X 

X 

) #F 

X 

- 21294621/468540250 SIN ( X X ) *E 9418697/ 


463 440'’40 

rps ( 

X 2 2 

X X ) »E * 1/468 X SIN ( 3 X ) ♦ 1/312 X COS C 

3 X 


) > 


LX 


F X 


^ 2 

X/ ? X ♦ X - 1 


- 4/7 X 

X/2 X ♦ ?/5 SIN 1 4/7 X J #E 


ORIGINAL PAGE IS 
POOR QUAUTV 


- 4/7 X 


9 r l/-» X ' 71402/459281 SIS C 4/7 X ) *E ♦ »353 74/ ?‘’<J640 S CO'i 


- 4/x X 

( 4/7 X ) fP 


S s 2 


- 4/7 X 

T = CHAIN ( X/2 X, - 73402/459231 SIN ( 4/7 X > *E *■ 2353X6/ 


- 4/ 7 X 

2296404 COS ( 4/7 X ) *« ) 
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2 

LX»5X +3X-1 


*» 2 X 

FX = X STN { 3 X ) »E 


2 2 X 2 2 X 

R = - 20/f»l51 X SIN < 3 X ) »£ - 6<3/5161 X CQS ( 3 X > »E 


2 X 

4- X6620S/2652592 1 X SIN ( 3 X ) «E ♦ 134700/26635921 X COS t 3 X ) 


? X 2 X 

- 659109350/1374679SS2S1 SIN { 3 X ) *E ♦ 441581922/137 


2 X 

457058281 COS ( 3 X ) *F 


5 = 1 

2 2 X ? 2 X 

T = - 20/5161 X SIN ( 3 X } 0E - 69/5161 X CQS C 3 x ) »£ 


2 X 

♦ 356205/26S35921 X SIN ( 3 X ) ♦ 134700/26635921 X CQS ( 3 X ) 


2 X 2 X 

• = - 650109350/137467983281 SIN ( 3 X ) *E «■ 441501922/137 


2 X 

467088281 COS ( 3 X ) 


ENO QF Q.O.F. SOLVER 
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ANALYSIS OF ALGORITHM 

Following the discussion in [1], we assume that there are functions M(n) and 
D(n), where M(n) - time to multiply two nth degree polynomials and D(n) - time 
to divide a polynomial of degree at most 2n by a polynomial of degree n, time being 
measured in terms of the number of bit operations required on a computer. In 
addition, we assume that the running time of the program is essentially deter- 
mined by the time to perform operations of the above types, so that we shall ignore 
the time-requirements of scalar multiplication and symbolic differentiation in our 
analysis. From this assumption it follows that the running time is determined by 
step 2 of the algorithm. Consequently, we restrict our attention to that step . 

We further restrict ourattention to case 2, where f(x) ■ x** exp(ort-i6)x, 3 0. 

2 2 2 

To facilitate the analysis, we set N(t) ■ t -2at+ot +6 , ■ deg L(t), and 

assume for convenience that iU2r ^ 3, where r is the multiplicity of oiilB 
as a root of L(t). We assume further that the method in Proposition 2 is used 
to implement step 2. In addition, we set q^ - deg B^, where AqCc), BQ<t) are 
polynomials as described in Proposition 2, with L » L/N*^. 

Theorem The time to execute step 2 of the algorithm in case 2 is bounded 
above by 

(r+2) D[^?] +D[1] +M[A-2r-2] + 21og2(2n+l) M[ (2n+l) (qQ+2)-(«-2r) ] , 
where Y denotes the least integer > Y. 


Proof Step 2 may be divided into three parts: 


a. 

determine 

L(t) - 

L(t)/N(t)’^, where N(t)’^|L(t), N(t)*"^^ | L(t) ; 

b. 

calculate 

Ao(t), 

Bo(t) 

such that AqL + BqN ■ I; 

c. 

calculate 

A^(t). 

B^(t) 

(2®1 

recursively by A . , ■ A [I*^B N' '], 

X 8 S 


2 

B 

s+l s 

until 

2® > 

n+l. 


We proceed to determine the running time for each step. 
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2 

a. This part involves computing successively = L/N, L 2 = Lj^/N ■ L/N , 

L = L ,/N - L/N*^, and terminating with L r N, which has a non-zero 
r I “ i, r 

remainder. The time to compute is bounded above by Since 

deg > deg > deg > . . . , the times to execute the remaining steps are 
all bounded above by D[is£,] as well so that the running time for part a is 
bounded by (r+1) . 

b. Since 4-2r ^ 3 by assunq)tion* we have that deg L > 3, hence L - N + T^, 

for polynomials S^, Tq, with deg T^ < 1. If deg T^ - 1, then N - S^^Tq + T^^ 
for polynomials T^, with deg Tj^ « 0, necessarily non-zero. It follows 

that AqL + B^N - I, for - "S^/T^, . ^S^Sq+I)/T^. Now 

deg Sq ■ 2,-2r-2 ^ 1 and deg ■ 1. It then follows that the work to calculate 
Aq and Bq In this situation is bounded above by D[* 5 S,-r] + D[l] + M[£-2r-2]. 

If deg Tq ■ 0, we may similarly obtain upper bounds on the work to calculate 
and on the degrees of Aq and Bq. in all cases the bounds do not exceed the 
corresponding ones obtained above. 


c. Set p = deg A , q » deg B . It follows readily from the recursive definitions 
s s s s 

c c 

of A and B that q “2 . p - q + 2 - (£-2r) » 2“(q^+2)- (H-2r). Hence 

S S S V s s u 

Pq £ 1 ?2 5 . • • • £ Pg ‘Jq — ^1 — ^2 — ' * * — perform 

each step of the recursion is bounded above by the time for the last step. 

Suppose now that s is the smallest integer such that 2® ^ n+1. Then the 

time to compute A is determined by the time to compute A , B n'’ To 

s ^ s-1 s-1 

get upper bounds on the degrees of these polynomials, note first that 

n+1 £ 2® £ 2n+l, by minimality of s. It then follows that p , £ *5(2n+l) (q_+2) 

a-1 ^ 

- (fc-2r), £ ^(2n+l)qQ, deg N' £ 2n+l. Adding the bounds, it follows 

that Pg_i *lg_i ^ £ (2n+l)(qQ+2) - (2,-2r). Thus the time to compute 

\-l \-l ^ ^® bounded above by 2M[ (2n+l) (qQ+2) - (i-2r)], and the total 
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time to compute 1. bounded ebove by 2. M[(2n+1) (,^+2) - (2-2r)! 

< 2 log 2 ( 2 n+l)M[( 2 n+l)(qQ+ 2 ) - (i-2r)J, since 2® < 2n+l. 

Having obtained bounds for the execution times of parts a, b» and 
step 2 , we obtain the overall bound announced In the theorem by adding 
together. Q.E.D. 


c of 
these 
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COMPARISON WITH OTHER METH(R)S 


Linear differential equations of the type considered here are commonly 
solved by the method of undetermined coefficients or the method of Laplace trans- 
forms. The former method requires In general r+n divisions of the polynomial 
L(t) over the complex n\imbers> while our method requires r+2 divisions over 
the reals. The latter method requires the factorization of the polynomial L(t) 
over the rationals, and In general will not give explicit answers unless all 
factors have degree £ A. See [5] and [A] for details on the two methods. 

Both methods have been Implemented as programs In the MACSYMA system, the 
former by Ivle [7] and the latter by Bogen [3]. It would clearly be of interest 
to compare the methods with ours in a practical sense, by Implementing the three 
in a common system and trying them on a common set of problems. 
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